#ifndef am_TempInterp_H
#define am_TempInterp_H

#include <string>
#include <tbb\concurrent_vector.h>
#include <boost\shared_ptr.hpp>

#include <borealis\Raster.hpp>
#include <borealis\AsciiRaster.hpp>
#include <borealis\Interpolation.hpp>
#include "..\AuroraTypes.hpp"

namespace Aurora
{
	/*
	Class: am_TempInterp
		Interpolates temperature over a spatial domain. 
	
	In Variables:
		Raster<float>	DEM variable
		string			Temperature header value from the met files
		string			Date header value from the met files
		string			Interpolation method
							- IDW => Inverse Distance Weighting
							- TPST => Thin plate spline with tension = 0.1
		string			Lapse Rate method
							- LLRA_const => Linear Lapse Rate : constant value of 6.5oC/km
							- LLRA_var   => Linear Lapse Rate : variable value for each month following Liston  (2006)
							- NSA_const
							- NSA_var
	Out Variables:
		string			Name of gridded Temperature

	Reference:
		Liston, Glen E., and Kelly Elder. 2006. A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). Journal of hydrometeorology 7: 217-234.

	*/
	class am_TempInterp : public AuroraTypes::ModuleBase
	{

	
		class RaiseInterpTemps: public Borealis::RasterVisitor<float>
		{
			float m_lapse_rate;
			int m_month;
			std::string m_interpMethod;
			boost::shared_ptr<Borealis::Raster<float> > m_dem;

		public:

			RaiseInterpTemps(const float lapse_rate,std::string interpMethod,boost::shared_ptr<Borealis::Raster<float> > dem,int month )
				: m_lapse_rate(lapse_rate)
			{
				m_interpMethod = interpMethod;
				m_dem = dem;

			}

			void operator()(float* cell, unsigned int row, unsigned int col)
			{
				if(m_interpMethod == "LLRA_const")
				{
					*cell = *cell+m_lapse_rate*(0-m_dem->GetRasterAt(row,col));//this is because we correct to sea level = 0.0	so negate the lapse rate (we're going up)
				
				}
				else if(m_interpMethod == "LLRA_var")
				{
					switch(m_month)
					{
					case 1:
						m_lapse_rate=0.0044F;
						break;
					case 2:
						m_lapse_rate=0.0059F;
						break;
					case 3:
						m_lapse_rate=0.0071F;
						break;
					case 4:
						m_lapse_rate=0.0078F;
						break;
					case 5:
						m_lapse_rate=0.0081F;
						break;
					case 6:
						m_lapse_rate=0.0082F;
						break;
					case 7:
						m_lapse_rate=0.0081F;
						break;
					case 8:
						m_lapse_rate=0.0081F;
						break;
					case 9:
						m_lapse_rate=0.0077F;
						break;
					case 10:
						m_lapse_rate=0.0068F;
						break;
					case 11:
						m_lapse_rate=0.0055F;
						break;
					case 12:
						m_lapse_rate=0.0047F;
						break;
					//no default case here as we can't throw from within a parallel_for loop
					//if we get this far, we already know that we have the right months as it was processed below
					//which would have thrown an exception on an error. 
					}
					*cell = *cell+(m_lapse_rate)*(0-m_dem->GetRasterAt(row,col));//this is because we correct to sea level = 0.0	so negate the lapse rate (we're going up)	
				}
				else if (m_interpMethod == "NSA_const")
				{
					double Po = 100000.0; //sea level pressure (Pa)
					double Cp = 1005.0; //specific heat of dry air J/[KgK]
					double R = 8.3143; //gas constant J/(molK)
					double m = 0.02897; //molecular weight of dry air kg/mol
					double Tb = 300.0; //assumed sea level temperature 300K
					double elev = m_dem->GetRasterAt(row,col); //cell elevation
					double g = 9.810616; //gravity m/s/s
					double lapse_rate = 0.0065; //K/m  negative here, going up

					//pressure at cells' elevation
					double Pz = Po * pow((double)(Tb/(Tb+lapse_rate*elev)),(double)((m*g)/(lapse_rate*R)));
					
					double Theta = *cell;

					double ratio = (Po/Pz);
					double exp = R/(m*Cp);

					double Ta = ( Theta/pow(ratio,exp) );
					Ta -= 273.15;
					*cell = Ta;

				}
			}
		};
		
		

	public:

		//in variables:
		//DEM
		//temp
		//date
		//interp method (IDW is only supported)
		//lapse rate method ("LLRA_const", "LLRA_var")

		int DEM;
		int TEMP;
		int DATE;
		int INTERP;
		int LAPSERATE;

		//outvariable
		//Interpolated Temperature raster
		int FILE;

		am_TempInterp(std::string n)
		{
			 DEM		=0;
			 TEMP		=1;
			 DATE		=2;
			 INTERP		=3;
			 LAPSERATE	=4;

			 FILE		=0;

			 name = n;
		}

	

		void run(AuroraTypes::StationList* stations,AuroraTypes::VariableList* variables, AuroraTypes::InVariables& inVars, AuroraTypes::OutVariables& outVars)
		{

			float lapse_rate=0.0065F;
			std::string m_interpMethod;
			boost::shared_ptr<Borealis::Raster<float> > m_dem;
			boost::shared_ptr<Borealis::AsciiRaster> m_temperature ;


			AuroraTypes::VariableList::const_accessor a;

			//bail if we don't have any met stations
			if(stations->size() == 0)
				throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) +"No met data");
			
			try
			{
				m_dem = GetRaster(variables,inVars[DEM]);

			}
			catch (std::runtime_error e)
			{
				 throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) +" - failed to get DEM parameter");
			}
			

			m_temperature.reset(new Borealis::AsciiRaster(m_dem->GetRow(),m_dem->GetCol()));
			m_temperature->SetXllCorner(m_dem->GetXllCorner());
			m_temperature->SetYllCorner(m_dem->GetYllCorner());
			m_temperature->SetCellSize(m_dem->GetCellSize());

			m_interpMethod = inVars[INTERP];
			//for each station lower the temperature to a common level z
			//build a new point and save it for the interpolation scheme
			tbb::concurrent_vector<Borealis::RasterPoint<float> > points;
			std::string lapseRate = inVars[LAPSERATE];


			//get the date from the first station
			//we assume at this point that all the stations are the same length, have the same time step, etc
			int month = (stations->at(0))->Now().Month();

			//correct the reference temperatures
			for(unsigned int i = 0; i < stations->size();i++)
			{
				double temperature = stations->at(i)->Now().Get<float>(inVars[TEMP]);
				unsigned int col = stations->at(i)->GetX();
				unsigned int row = stations->at(i)->GetY();
				
				float newTempterature=-9999.0;

				//constant lapse rate
				if(lapseRate == "LLRA_const")
				{
					lapse_rate=0.0065F;		
					newTempterature =  temperature - lapse_rate*(0.0 - stations->at(i)->GetElevation());
				}
				else if(lapseRate == "LLRA_var")
				{
					switch(month)
					{
					case 1:
						lapse_rate=0.0044F;
						break;
					case 2:
						lapse_rate=0.0059F;
						break;
					case 3:
						lapse_rate=0.0071F;
						break;
					case 4:
						lapse_rate=0.0078F;
						break;
					case 5:
						lapse_rate=0.0081F;
						break;
					case 6:
						lapse_rate=0.0082F;
						break;
					case 7:
						lapse_rate=0.0081F;
						break;
					case 8:
						lapse_rate=0.0081F;
						break;
					case 9:
						lapse_rate=0.0077F;
						break;
					case 10:
						lapse_rate=0.0068F;
						break;
					case 11:
						lapse_rate=0.0055F;
						break;
					case 12:
						lapse_rate=0.0047F;
						break;
					default:
						throw std::runtime_error("In " + std::string(__FILE__) +" @ line " + boost::lexical_cast<std::string>(__LINE__) + "unknown month for lapse rate calculation");				
					}
					newTempterature =  temperature - lapse_rate*(0.0 - stations->at(i)->GetElevation());
				}
				else if(lapseRate == "NSA_const")
				{
					double Po = 100000.0; //sea level pressure (Pa)
					double Cp = 1005.0; //specific heat of dry air J/[KgK]
					double R = 8.3143; //gas constant J/(molK)
					double m = 0.02897; //molecular weight of dry air kg/mol
					double Tb = 300.0; //assumed sea level temperature 300K
					double elev = stations->at(i)->GetElevation(); //station elevation
					double g = 9.810616; //gravity m/s/s
					double lapse = -0.0065F; //K/m

					//pressure at station's elevation
					double Pz = Po * pow(Tb/(Tb+lapse*elev),(m*g)/(lapse*R));
					temperature +=273.15; //to K
					double ratio = (Po/Pz);
					double exp = R/(m*Cp);
					double theta = temperature * pow(ratio,exp);
					newTempterature =  theta;
					lapse_rate = lapse;
				}

				Borealis::RasterPoint<float> point(row,col,newTempterature);
				points.push_back( point );
			}

			if(m_interpMethod == "IDW")
			{
				RaiseInterpTemps* visitor =new RaiseInterpTemps(lapse_rate,lapseRate,m_dem, month);
				Borealis::Interpolation::InverseDistanceWeighting(points, m_temperature, visitor);
				delete visitor;
			}
			else if(m_interpMethod == "TPST")
			{
				RaiseInterpTemps* visitor =new RaiseInterpTemps(lapse_rate,lapseRate,m_dem, month);
				Borealis::Interpolation::ThinPlateSplineT(points, m_temperature, visitor);
				delete visitor;
			}
			else
			{
				throw std::runtime_error("Unknown interpolation method");
			}
			
			AuroraTypes::VariableList::accessor a2;
			variables->insert(a2,outVars[FILE]);		
			a2->second = m_temperature;

		}

	};
}

#endif